%MY_LAPLACEEDGE Summary of this function goes here
% This algorithm is divided into three steps
% Step 1 calculates two gaussians and takes the differences between them,
% these results in an image with reduced spatial information, keeping
% intact the greater constrast differences (edges)
%
% Step 2 smoothes the filtered image on step1 by zeroing out any high
% constrast pixel that is not really and edge. We do this by comparing the
% pixel with its neighbors and checking if it really is between a negative
% (background) and a positive (foreground) region
%
% Step 3 removes any noise that still remains through the pre-stipulated
% threshold
function [ Gthres, Gsmooth, Gmh  ] = my_marrhildrethEdge( img, sigma0, sigma1, threshold )

h = size(img,1);
w = size(img,2);

% Step 1 - difference of gaussians
Gmh = my_differenceOfGaussians(img, sigma0, sigma1);
Gsmooth = zeros(h,w);
Gthres = zeros(h,w);

% Step 2 - zero crossings
for i = 2:h
    for j = 2:w
        if Gmh(i,j) > 0
            if (Gmh(i,j+1)>=0 && Gmh(i,j-1)<0) || (Gmh(i,j+1)<0 && Gmh(i,j-1)>=0)
                Gsmooth(i,j) = Gmh(i,j+1);
            elseif (Gmh(i+1,j)>=0 && Gmh(i-1,j)<0) || (Gmh(i+1,j)<0 && Gmh(i-1,j)>=0)
                Gsmooth(i,j) = Gmh(i,j+1);
            elseif (Gmh(i+1,j+1)>=0 && Gmh(i-1,j-1)<0) || (Gmh(i+1,j+1)<0 && Gmh(i-1,j-1)>=0)
                Gsmooth(i,j) = Gmh(i,j+1);
            elseif (Gmh(i-1,j+1)>=0 && Gmh(i+1,j-1)<0) || (Gmh(i-1,j+1)<0 && Gmh(i+1,j-1)>=0)
                Gsmooth(i,j) = Gmh(i,j+1);
            end
        end
    end
end

% Step 3 - Threshold to remove noise
Gopt = im2uint8(Gsmooth);
Gthres = Gopt > threshold;


end

